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Abstract 

In this paper, we describe a numerical continuation method that enables har- 
monic analysis of nonlinear periodic oscillators. This method is formulated as 
a boundary value problem that can be readily implemented by resorting to a 
standard continuation package - without modification - such as AUTO, which 
we used. Our technique works for any kind of oscillator, including electronic, 
mechanical and biochemical systems. We provide two case studies. The first 
study concerns itself with the autonomous electronic oscillator known as the 
Colpitts oscillator, and the second one with a nonlinear damped oscillator, 
a non-autonomous mechanical oscillator. As shown in the case studies, the 
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proposed technique can aid both the analysis and the design of the oscilla- 
tors, by following curves for which a certain constraint, related to harmonic 
analysis, is fulfilled. 

Keywords: continuation methods; harmonic analysis; nonlinear oscillator; 
boundary value problem. 



1 Introduction 

Continuation methods provide an efficient tool for analyzing systems of non- 
linear algebraic equations whose solutions form a one-dimensional continuum. 
When dealing with periodic solutions of systems of ordinary differential equa- 
tions (ODEs), we continue solutions by solving a boundary value problem 
(BVP). We can either solve this BVP explicitly, when possible, or implicitly, 
by using numerical continuation tools. 

The BVPs that are used to continue "standard" objects, such as equilib- 
rium points, periodic, homoclinic, and heteroclinic orbits and their bifurca- 
tions in dynamical systems are well-known (see, for instance, [Kuznetsov, 2004" 



for an overview). However, it is often advantageous to formulate new BVPs to 
continue "non-standard" objects, such as invariant manifolds [Doedel et ai, 2006 



slow manifolds [Desroches et ai, 2008| , and coherent structures such as spi- 
ral waves and other defects in oscillatory media [Bordyugov fc Engel, 200T 
Champneys fc Sandstede, 2007| . 



In a recent paper [Cochelin &: Vergez, 2009 , a combination of harmonic 



analysis and continuation techniques (based on the harmonic balance method) 
was proposed. The main result is a new numerical continuation tool that 
provides standard results of continuation analysis after a preliminary refor- 
mulation of the problem in terms of harmonic balance. 

In this paper, we show how to define a novel BVP that enables har- 
monic analysis by using standard numerical continuation tools. The pro- 
posed BVP allows, in the spirit of [Doedel et ai, 2006 Desroches et ai, 2008 



Bordyugov fc Engel, 20071 [Champneys ^ Sandstede, 2007|, for "non-standard" 



continuations focused on selected harmonic components of a solution without 
an ad hoc simulator like the one described in [Cochelin &: Vergez, 2009 . 



The first application in this paper deals with the design of an autonomous 
electronic circuit (the Colpitts oscillator). Its aim is to obtain parameter 
charts that can be immediately understood by designers of electronic oscil- 
lators. In particular, these charts show how the limit cycle corresponding to 
a periodic regime changes in a familiar circuit parameter plane. The guide- 
lines for the designer are provided, for instance, by curves corresponding to 
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solutions with fixed amplitudes of some harmonic components or with a fixed 
ratio of two harmonic components. 

The second application entails the harmonic analysis of a non-autonomous 
mechanical oscillator (nonlinear damped oscillator). In particular, we use the 
proposed BVP to analyze the so-called jump phenomenon [Schmidt fc Tondl, 1986 . 

The main advantages of the proposed approach are the following: 

• the BVP is formulated in such a way that the public-domain software 
package AUTO-07P [Doedel fc Oldeman, 2009] can solve it; 

• compared to simulations, computation times are generally lower, since 
numerical continuation packages operate directly on system invariant 
sets; 

• the proposed procedure is reasonably easy to use by those who are 
familiar with the analysis of nonlinear dynamical systems. 

On the other hand, one needs to take care of the following aspects: 

• for complex systems, such as a realistic radio-frequency electronic os- 
cillator, the procedure is effective only if preceded by a modeling phase 
where one defines a suitable model of the oscillator [Bizzarri et ai, 2009 
— such a model should be as simple as possible, but able to capture 
the essential features of the system; 

• those who are not familiar with numerical continuation packages require 
a preliminary training. 

This paper is organized as follows. Section |2] summarizes the basic el- 
ements that the proposed technique is based on. The BVP formulation is 
described in Sec. [31 whereas Sees. H] and O are devoted to two case studies. 
Some conclusions are drawn in Sec. El 

2 Basic elements 

Let the following system of ODEs describe an oscillator, which can be either 
autonomous or non-autonomous with periodic forcing: 

dx 

x = — = g(x, t;p), x,g e M", p G M", t G M. (1) 
at 

When this oscillator reaches a periodic regime, it produces signals of generic 
period T, that is, frequency / = ^ and angular frequency u = 2nf. The 
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Fourier series expansion of each state variable is 



Xj(t) = aoj + [ttkj sm{kut) + bkj cos{kut)] , (2) 



k=l 



where the index j G {l,...,n} selects the state variable, a^j is the mean 
value of Xj(t) over T and we obtain the other coefficients by projecting Xj(t) 
on the corresponding basis functions 



2 /•*+^ 

ttkj = — / Xji^r) sm{kujT)dT 

(3) 



2 



t+T 



bkj = — / Xj{T) cos{kuiT)dT 
T Jt 

In general, continuation methods show how system invariants (for ex- 
ample, equilibrium points or limit cycles) depend on one or more control 
parameters. One of the key elements of continuation theory is that invariant 
sets and their bifurcations are revealed when so-called test functions equal 
zero. Many test functions are included in the most diffused numerical con- 
tinuation packages [Doedel fc Oldeman, 2009 Dhooge et ai, 2003 , but, of 



course, user-defined test functions can be added. For instance, in this paper 
we define test functions of the kind /(S'a, S*;,, T, /Tref); where Sa and St de- 
note a subset of and Ub Fourier coefficients {atj} and {bkj}, respectively, 
and -R'ref is a constant reference value. 

We can use this formulation to solve different kinds of problems. Firstly 
we can continue the amplitude of a harmonic component with respect to a 
single parameter. However, of further interest are continuations with respect 
to two parameters, adding a further constraint. For instance, iso-harmonic, 
iso-ratio, and iso-energy continuations are possible. For iso-harmonic contin- 
uations we fix the amplitude of a harmonic component, for iso-ratio continua- 
tions the ratio between amplitudes of different harmonics, and for iso-energy 
continuations a sum of squares of a limited number of harmonic amplitudes 
representing almost the whole power spectrum of the analyzed signal. 

In the next section, we set up the EVP that is solved for these continua- 
tions. 



3 Setup 

We start by deriving a simplified model of the oscillator given by a small sys- 
tem of ODEs, as is usual when dealing with both the analysis and synthesis 
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of a dynamical system. We call this system the original system. By assuming 
some reasonable modeling hypotheses we obtain these equations, which are 
expressed in terms of state variables. It is advantageous, but not compulsory, 
to normalize these equations and shift the origin of the normalized state space 
to a "significant" equilibrium point. Such an equilibrium is stable for some 
parameter configuration, and, by varying one parameter, undergoes a super- 
critical Andronov-Hopf bifurcation, which marks the appearance of a family 
of asymptotically stable periodic solutions evolving around the (unstable) 
equilibrium. 



3.1 The original system 

For autonomous oscillators, the original system is simply the ODE system 
modeling the oscillator, given by i; = g{x(t);p). For a non-autonomous 
oscillator with periodic forcing, we can obtain an equivalent autonomous 
oscillator by adding a nonlinear oscillator with the desired periodic forcing 
as one of the solution components (see, for instance, the AUTO-07P demo 
frc [Doedel fc Oldeman, 2009 Alexander et al., 1990| ). In particular, for a 
sinusoidal forcing we can use the following secondary oscillator (very close to 
the normal form of the supercritical Andronov-Hopf bifurcation): 



V = av + f3w + f (f ^ + w'^) 
w = —pv + aw + w[v + w ), 

which for a < asymptotically converges to the origin and for a > has the 
asymptotically stable solution v = sin(/3t), w = cos(/3t). For instance, if the 
first state variable of the original system obeys the differential equation Xi = 
gi{x{t)]p) + ccos{ut) and the state of the original system is x = [xi, . . . 
then in the corresponding autonomous system the equation for the first state 
variable is given by Xi = gi{x{t)]p) + cw, where [5 = u and the state vector 
is redefined as a; = [xi, . . . , a;„, ty]. Hence, if the original system is a non- 
autonomous oscillator with periodic forcing, we can also recast it to the 
autonomous system x = g{x{t);p). 

We can analyze equilibria, their Hopf bifurcations and emanating limit 
cycles in the possibly recast original system. The original system, together 
with some information from the limit cycle close to the Hopf bifurcation, is 
then used to construct and initialize the full continuation system defining the 
EVP. 
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3.2 Initialization 



As an initial solution for the BVP problem, we could take a pre-computed 
periodic orbit, carry out a signal analysis to find the coefficients of its Fourier 
expansion, and substitute that into the system. Another way (which we 
describe here) is to use standard continuation techniques provided by AUTO, 
following an equilibrium point that undergoes an Andronov-Hopf bifurcation. 

In the non-autonomous case with sinusoidal forcing, we can easily find 
this bifurcation by varying the parameter a in Eqs. (jl]) across zero. 

The limit cycle that emanates from the Andronov-Hopf bifurcation can 
then be continued a small distance away from the bifurcation, where it will 
still have the approximate form 



Here u is the purely imaginary part of the corresponding eigenvalue of the 
equilibrium, which we can obtain using the period T that AUTO provides: 
bj = 211 /T. Now the Fourier coefficients can be trivially derived, comparing 
Eqs. ([2]) and ([5]) : 



3.3 The full continuation system 

The full continuation system is given by the following equations: 

• Non-autonomous differential equations: 

x = Tg{x{t),t-p), x,^?gM", peW, t,TeR 

i = l ^ ' 

Here the possibly recast original system Eq. ([1]) is rescaled so that 
a T-periodic solution of Eq. ([T]) is a 1-periodic solution for the first 
equation in Eqs. ([7]). The non-periodic equation i = 1 is added to 
make the system look autonomous to continuation software. 

• Boundary conditions that define a periodic orbit on the t-interval [0, 1]: 



x{t) = A sin ut + B cos ut. 



(5) 



aoj — 0, 

aij = A = Xj{0)/uj, 
a-kj = 0, 



hj = B = a;,(0), 

bkj = 0, k> 1, j = l,...,n 



(6) 



a^i(O) = j = l,...,n 

t(0) = 
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• Integral conditions: 



1 

x{t)xo\d{t)dt = 



/ {2xj{t) sin(27rH) — akj)dt = for any atj E Sa 
Jo 

/ (2xj(t) cos(27rA;t) — hkj)dt = for any bkj E Si, 
Jo 

f f{Sa^St,T,KnEF)dt = 

Jo 



(9) 



Here the first condition is the standard integral phase condition [Kuznetsov, 2004" 



where Xoid(^) denotes the previous point on a continuation branch, 
the second and third conditions compute or fix the Fourier coefficients 
{ttkj} G Sa and {bkj} € Sb, and the fourth condition computes or fixes 
-^REF given the Fourier coefficients in Sa and Sb- 

This gives us a system of n + 1 ordinary differential equations, n + 1 
boundary conditions and Ua + Ub + 2 integral conditions. Adding a standard 
pseudo-arclength condition, this needs to be offset hj Ua + rib + S continuation 
parameters. 

A basic choice for the parameters is to continue a periodic orbit in one of 
the system parameters, the period T, and the Ua + rib + l values in Sa^ SbU 
{Khef}- Note that in this case these last Ua + Ub + 1 values are effectively 
measured through integral conditions, where explicit test functions such as 

akj= [ {2xj{t) sm{2TT kt))dt (10) 
Jo 

would suffice. 

However, the integral conditions become more powerful if we like to keep 
something fixed: for instance, by fixing -/^ref and freeing up one more sys- 
tem parameter we can continue iso-amplitude curves. Moreover, in existing 
continuation software it is easier and arguably more elegant, if a little more 
computationally expensive, to stick to one full system with the same inte- 
gral conditions, than to switch between test functions and equivalent integral 
conditions. 
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4 Case study 1: an autonomous electronic os- 
cillator 



Increasing demand, in modern communication systems, of high-performance 
low-power radio-frequency circuits is driving the development of accurate 
simulation tools at the design stage. The simulation is particularly critical 
in the case of analog circuits (free- running or voltage-controlled oscillators), 
because of time consumption, and in the case of mixed-signal analog-digital 
circuits (frequency synthesizers or phase-locked loops), due to the modelling 
difficulties. In all cases, the analysis is usually non-trivial and performance 
verification requires extensive simulation. This is even more important when 
the main goal is to find how changes in circuit parameters (for example, 
the amplitude or frequency of an input generator, or a linear-element value) 
affect circuit performance. 

The main research lines in this field are twofold. Firstly there is the de- 
velopment of algorithms that speed up and make circuit simulations more 
reliable [Brambilla et al., 2005 Brambilla fc Storti-Gajani, 2008| . Secondly 



there is the study of methods that, starting from a simplified version of the de- 
signed circuit, aim to provide design criteria avoiding brute- force simulations 
of complex integrated circuits. Most of these methods are based on well- 
established theories such as harmonic balance, Volterra series (for weakly 
nonlinear oscillators), and bifurcation analysis of nonlinear dynamical sys- 
tems. Harmonic balance is used for an accurate determination of nonlinear- 
circuit operation bands Rizzoli fc Neri, 2009 Suarez et ai, 2006 and, asso- 



ciated with other methods, for calculating the periodic response of nonlinear 
dynamical circuits [Buonomo fc Lo Schiavo, 2003| . Volterra series are used 
for the analysis of nearly sinusoidal nonlinear oscillators |Hu et al., 1989 



Huang &: Chu, 1994| . Bifurcation analysis is used to optimize the design 



of oscillators [Maggio et al., 1999| and the locking range of injection-locked 
frequency dividers [Ghahramani et al., 2007| . A combination of harmonic- 
balance simulators and bifurcation control is exploited to obtain bifurcation 
control in microwave circuits, thus presetting the operation bands of com- 
plex circuits, such as synchronized and voltage-controlled oscillators and fre- 
quency dividers [CoUado Suarez, 2005| . This combination is also used to 
obtain robust and efficient oscillator analysis techniques (see, for example, 
[Bonani fc Gilli, 1999 Genesio et al, 1993 Gourary et al, 2000| and refer- 



ences therein). 

[Figure 1 about here.] 
Our case study is the Colpitts oscillator, shown in Fig. [H whose dynamics 
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were analyzed in detail in [Maggio et ai, 1999 De Feo et ai, 2000 De Feo &: Maggio, 2003 
The following simplified model, that we can easily substitute into the system 
defined in Section [3l adequately describes this oscillator: 



G 



X 



y 



{e-y-l)+z] 



(11) 



Q7(l-7) , ^ N 1 



where the system parameters are initialized as follows: 



Q 



R 



0.8; 7 



G 



G1 + G2' " VTRiG, + G2) 
Gi = G2 = 10-^ [F]; L = 10~'^[H]. 



2: 



(12) 



Moreover, we assume a^? = 1, that is, we assume the CB short-circuit 
forward current gain of the BJT transistor to be ideal. The variable z denotes 
the inductor current normalized with respect to Iq. The variables x and y 
are the voltages across Ci and G2 normalized with respect to = 25.9 mV 
(that is, the thermal voltage at room temperature). Time is normalized with 
respect to Tq = ^/LC^C^jtc^TCi) . 

Continuing the equilibrium at as the parameter G goes downwards from 
2 we find a Hopf bifurcation at G = 1. For G slightly greater than 1 there 
then exists a limit cycle for which the Fourier coefficients ai2, 612, and -R'ref 
are given by (see Eqs. (jHD) 



a^2 = m/C^^/T) 
&i2 = y(0). 



Gz{Q) 



/(27r/T), 



(13) 



12 



/)2 
"12- 



We focus on the state variable y since it corresponds to the voltage used as 
output of the oscillator. Moreover, we focus on its first harmonic component 
since we want the generated oscillation to be nearly sinusoidal. Thus we need 
to verify that higher order harmonics have small amplitudes when compared 
with the first harmonic component. 

We then follow the periodic orbits as a boundary value problem in (G, T, 012, 
6i2,i^^REF) with user-defined labels for the following values of K-^-ef'- 1? 4, 7, 
10, 13, 16, 19, 22, 25. Using these values as starting points we can find 
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iso-harmonic curves in [G, 7)-parameter plane for fixed values of -RTref by 
continuing in (Q, G, T, a^, b^). 



[Figure 2 about here.] 

We carried out the continuation in the parameter space {Q, G, T, au, 612), 
but show the results in Fig. [2] on the plane {R,Iq), according to the defini- 
tions of Q and G, to simplify the interpretation of the results from a circuit 
point of view. We excluded the parameter space region characterized by 
the presence of complex dynamics [Maggio et al., 1999 De Feo et al., 2000] 



from our analysis and focused on the region characterized by the presence of 
a stable "simple" harmonic cycle. 

In the upper panel, along the black dashed curves the period (normalized 
with respect to Tq) is constant: 6.3 for the lower curve, and 8.7 for the 
upper curve, with a step of 0.3. In the lower-left panel, along the black 
solid curves the amplitude Ai (normalized with respect to Vr) of the first 
harmonic component of y is constant (1 for the lower curve, 25 for the upper 
curve, step of 3), whereas along the black dashed curves the amplitude A2 
(normalized with respect to Vr) of the second harmonic component of y is 
constant (1, 5, 10, and 15 from the lower curve to the upper one). 

In the left-right panel, we follow the periodic orbits as a boundary value 
problem in (G, T, ai2, 612, ^22, ^22, -^^ref)- Along the black solid curves the 
ratio ^ (= K^ef) is constant: 3 for the upper curve, 5 for the lower curve, 
with a step of 1. The interpretation of these results is straightforward: for 
instance, if we want to fix both circuit parameters to have an oscillation fre- 
quency / = To/6.3, and ensure low values of R to increase the quality factor 
Q, we can just properly adjust Iq along the lower curve in the upper panel of 
Fig. [2j The results shown in the lower panels provide useful information if a 
periodic signal whose first harmonic component is predominant with respect 
to the second one is of interest. 

Of course, we can obtain other iso-curves besides those shown in Fig. [2J 
For example, we may chose values of Iq and R ensuring a desired ratio ^ 
for given p and q. 

Figure [3] shows the obtained iso-harmonic (with constant Ai) curves in 
the parameter subspace {R, Jq, T). Of course, the projection of the curves on 
the plane {R, Iq) gives the black sohd curves displayed in the lower-left panel 
of Fig. 121 

[Figure 3 about here.] 
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5 Case study 2: a non- autonomous mechanical 
oscillator 



The nonlinear damped oscillator is a model widely used to represent shock 
absorbers [Lang et ai, 2006 and is given by the following system (with di- 



mensionless variables and parameters): 

(14) 



X = — cos(ci;t) (cix + C2X^ + c^x^ + ky) 

m m 



y = x 

where m = 240, Ci = 296, Ca = 3000, cg = 800, and k = 240(47r)2. The 
reference angular frequency of the system is cuq = ^yk/m = An. 

We can easily recast the non-autonomous system as an autonomous sys- 
tem by adding the auxiliary oscillator (jlj) to Eqs. (fT4|) . For a < 0, the 
resulting system has a stable equilibrium at the origin, which undergoes a 
Hopf bifurcation for a = 0. For small positive a there then exists a limit 
cycle for which the Fourier coefficients ai2, and K^^f are approximately 
given by, and determined numerically as (see Eqs. (^) 

au = 2/(0)/(27r/T) = x(0)/(27r/T), 

^12 = 2/(0), (15) 

Kref = ^a?2 + ^12- 

We focus on the state variable y since it represents the position of the 
mechanical oscillator. 

Once the initialization is completed, we work with the full continuation 
system and continue the limit cycle with respect to (a, T, an, 611, Kref) until 
we reach a = 1, corresponding to the correct sinusoidal forcing. 

In this case, we monitor the amplitude of the first harmonic with re- 
spect to one of the parameters u and A (properly normalized), besides 
{T, ail, bii, Kref), to provide evidence for the presence of a "jump phe- 
nomenon". The "jump phenomenon" is a characteristic feature of many non- 
linear oscillators, where the response amplitude changes suddenly at some 
critical value of the excitation frequency [Schmidt fc Tondl, 1986| . Many 
widely used methods such as the Harmonic Balance (HBM) and Nonlinear 
Output Frequency Response Function (NOFRF) [Peng et ai, 2008] barely 
capture this phenomenon. Figure H] shows the amplitudes of the first three 
harmonics as functions of the normalized frequency u/uo (upper panels) and 
of the normalized amplitude A/m (lower panels). The results perfectly match 
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the reference diagrams reported in [Peng et ai, 2008 . 

[Figure 4 about here.] 

We remark that the proposed continuation method does not require any 
approximation, thus providing excellent accuracy. Moreover, it is based on 
tools (such as AUTO) that are reliable and widely tested. 

Finally, by performing continuations similar to those shown for the Col- 
pitts oscillator, we can also obtain simulations pertaining to the design of 
mechanical oscillators, that is, we obtain the values of parameters that en- 
sure a desired behavior. 

6 Conclusions 

We proposed a technique, based on numerical continuation, that enables 
harmonic analysis of a given nonlinear oscillator without resorting to any ap- 
proximation. Moreover, a designer can choose some oscillator parameters to 
obtain a desired behavior, by analysing some of the curves that are obtained 
by this technique. More realistically, since the model is an approximation of 
the real system, the proposed method provides at least reference values of 
the bifurcation parameters. More accurate simulations focused on restricted 
portions of the parameter space can then refine these values. 

The main advantage of this technique is that it enables the analysis of 
even relatively complex oscillators (both autonomous and forced) by using 
software tools that are reliable and optimized. This makes the development of 
ad hoc software for this kind of analysis unnecessary and ensures an excellent 
accuracy of the results. 

The main limit of this technique is that it requires a thorough knowledge 
of continuation methods and/or software packages for numerical continua- 
tion. Moreover, for oscillators forced by non-sinusoidal periodic signals, it 
can be non-trivial to define the auxiliary equations needed to make the sys- 
tem autonomous. In these cases, it may be more convenient to use the results 
of a signal analysis carried out on a pre-computed periodic orbit. The coef- 
ficients of the thus obtained Fourier expansion can be substituted into the 
system, which can then be used as an initial solution for the BVP problem. 
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Figure Captions 

1 The considered Colpitts oscillator. Vcc = 3 V, Ci = C2 = 1 
/iF, L = 10 niH 

2 Charts for the Colpitts oscillator on the parameter plane [R, /q) . 
The grey line is the Hopf bifurcation curve in all panels. Iso- 
period curves (upper panel), iso-Ai (solid) and iso-yl2 (dashed) 
curves (lower-left panel), and iso-ratio ^ (lower-right panel). . 

3 Iso-harmonic curves in the parameter subspace (i?, Jo, T). . . . 

4 Amplitudes of the first (first column), second (second column) 
and third (third column) harmonics with respect to u/uq (up- 
per panels) and to A/m (lower panels) 
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Figure 1: The considered Colpitts oscillator. Vcc = 3 V, Ci = C2 = 1 /iF, 
L = 10 mH. 
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Figure 2: Charts for the Colpitts oscillator on the parameter plane (i?, /q). 
The grey line is the Hopf bifurcation curve in all panels. Iso-period curves 
(upper panel), iso-Ai (solid) and iso-/l2 (dashed) curves (lower-left panel), 
and iso-ratio ^ (lower-right panel). 



18 




Figure 3: Iso-harmonic curves in the parameter subspace {R,Io,T). 
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Figure 4: Amplitudes of the first (first column), second (second column) and 
third (third column) harmonics with respect to oj/ujq (upper panels) and to 
A/m (lower panels). 
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